#Joystick model for tanks according to Bakalis et al. (2017)
 
import openseespy.opensees as op
import vfo.vfo as vfo
import os
import numpy as np
import matplotlib.pyplot as plt
import time
import Units as ut
import scipy as scp


op.wipe();

direc = './Results'

if(os.path.exists(direc)==False):
    os.mkdir(direc);
    
    
#Obtain ground motion parameters     

DT = np.loadtxt('./475/GMRs/GMR_dts.txt')
DUR = np.loadtxt('./475/GMRs/GMR_durs.txt')
#print(Npts)

for GM in range(0,len(DT)):

    GMfile = './475/GMRs/GMR_'+str(GM+1)+'.txt'
    dts = DT[GM]
    Npts = int(DUR[GM]/DT[GM])

    print('Model generation started...')
        
    R = 13.9*ut.m #Tank Radius in m
    H = 16.5*ut.m #Contained Liquid Height in m
    nsp = 8 #Number of beams used to modelling the base plate, Bakalis et al. (2017)
    ro = 998*ut.kg_m3 #liquid density in kg/m3 
    fillr = 0.95 #Pencentage of tank's height filled with liquid 

    m = ro*H*np.pi*R**2*fillr

    #Compute impulsive and convective masses, heights, and periods according to Calvi & Nascimbene (2023)
    n = np.linspace(0,245,245, dtype='int')
    gamma = H/R
    vn = np.zeros(len(n))
    for i in range(len(n)):
        vn[i] = np.pi*(2*n[i]+1)/2
        
    #Compute modified bessel function of the first kind and its derivative   
    I1 = np.zeros(len(n))
    I1p = np.zeros(len(n))

    for i in range(len(n)):
        I1[i] = scp.special.i1(vn[i]/gamma)
        I1p[i] = scp.special.i0(vn[i]/gamma)- scp.special.i1(vn[i]/gamma)/(vn[i]/gamma)

    #Impulsive mass
    mi = np.sum(2*m*gamma*(I1/(I1p*vn**3)))  
    #Impulsive height
    hi = H*np.sum((-1)**n*I1*(vn*(-1)**n-1)/(vn**4*I1p))/np.sum((I1/(I1p*vn**3)))
    #Coefficient for Impulsive period, depends on gamma, page 139 in Calvi and Nascimbene (2023)
    Ci = 6.2
    tw = 10.61*ut.mm  #wall thickness
    rom = 998   #liquid density in kg/m3, don't transform to ton/m3 becuase of formula for Ti
    Es = 210000000000 #Modulus of elasticity of steel, in Pa, don't transform to kPa becuase of formula for Ti 
    #Impulsive period
    Ti = Ci*H*(rom/Es)**0.5/(tw/R)**0.5
    #Equivalent impulsive stiffness
    ki = 4*np.pi**2*mi/Ti**2

    ksiI = 0.02 #Damping ratio of impulsive component 
    #Define damping constant of impulsive component
    ci = ksiI*2*(ki*mi)**0.5

    #Convective mass
    lam1 = 1.8412  #first root of the first derivative of the Bessel function of the first kind
    mc = m*2*np.tanh(lam1*gamma)/(lam1*gamma*(lam1**2-1)) 
    #Convctive height
    hc = H*(1+(1-np.cosh(lam1*gamma))/(lam1*gamma*np.sinh(lam1*gamma)))
    #Fundamental Convective period
    Tc = 2*np.pi*(R/ut.g)**0.5/(lam1*np.tanh(lam1*gamma))**0.5
    #Equivalent convective stiffness
    kc = 4*np.pi**2*mc/Tc**2

    ksiC = 0.005 #Damping ratio of convective component 
    #Define damping constant of convective component
    cc = ksiC*2*(kc*mc)**0.5

    op.model('basic', '-ndm', 3, '-ndf', 6) 

    #Define central node
    nodeTagCent = 1
    op.node(nodeTagCent,0,0,0)

    #Define impulsive mass nodes
    nodeTagi1 = 2
    nodeTagi2 = 3
    op.node(nodeTagi1, 0,0,hi)
    op.node(nodeTagi2, 0,0,hi)

    #op.equalDOF(nodeTagi1,nodeTagi2,3,4,5,6)

    #Define impulsive mass
    op.mass(nodeTagi2, mi,mi,0)

    #Define convective mass nodes
    nodeTagc1 = 4
    nodeTagc2 = 5
    op.node(nodeTagc1, 0,0,hc)
    op.node(nodeTagc2, 0,0,hc)

    #op.equalDOF(nodeTagc1,nodeTagc2,3,4,5,6)

    #Define convective mass
    op.mass(nodeTagc2, mc,mc,0)

    #Define nodes on tank circumference
    dth = 2*np.pi/nsp
    nodeTagCirc1 = 10
    nodeTagCirc2 = 100

    for i in range(0,nsp):
        op.node(nodeTagCirc1+i, round(R*np.cos(dth*i),5),round(R*np.sin(dth*i),5),0)
        op.node(nodeTagCirc2+i, round(R*np.cos(dth*i),5),round(R*np.sin(dth*i),5),0)
        op.fix(nodeTagCirc2+i, 1,1,1,1,1,1)
        #op.equalDOF(nodeTagCirc2+i,nodeTagCirc1+i,1,2,4,5,6)

    #Define Materials

    #Parameters of Pinching4 material model used to model the base plate uplift behaviour 
    matTag1 = 1

    ePf1 = 350
    ePd1 = 0.0038
    k1 = ePf1/ePd1

    ePf2 = 550
    ePd2 = 0.09
    k2 = (ePf2-ePf1)/(ePd2-ePd1)

    sigAct = ePf1
    beta =0.38
    epsSlip = 0.0
    epsBear = 0.2

    ePf3 = 1800
    ePd3 = 0.4

    k3 = (ePf3-ePf2)/(ePd3-ePd2)

    rBear = k3/k1

    #strain = [eNd4,eNd3,eNd2,eNd1,0,ePd1,ePd2,ePd3,ePd4]
    #stress = [eNf4,eNf3,eNf2,eNf1,0,ePf1,ePf2,ePf3,ePf4]
    #Define material of springs modelling the base plate uplift 
    op.uniaxialMaterial('SelfCentering',matTag1,k1,k2,sigAct,beta,epsSlip,epsBear,rBear)

    matTag2 = 2
    K = 800/0.3
    op.uniaxialMaterial('Elastic',matTag2,K,0,ut.Ubig)

    #Define impulsive and convective elastic materials
    ImpMatTagS = 3
    op.uniaxialMaterial('Elastic',ImpMatTagS,ki)
    ConvMatTagS = 4
    op.uniaxialMaterial('Elastic',ConvMatTagS,kc)

    #Define impulsive and convective viscous materials
    ImpMatTagD = 5
    op.uniaxialMaterial('Viscous',ImpMatTagD,ci,1)
    ConvMatTagD = 6
    op.uniaxialMaterial('Viscous',ConvMatTagD,cc,1)

    #Define viscoelastic materials for impulsive and convective components. 
    ImpMatTag = 7
    op.uniaxialMaterial('Parallel',ImpMatTag,ImpMatTagS,ImpMatTagD)
    ConvMatTag = 8
    op.uniaxialMaterial('Parallel',ConvMatTag,ConvMatTagS,ConvMatTagD)

    rigidMatTag = 9
    op.uniaxialMaterial('Elastic',rigidMatTag,ut.Ubig)

    flexMatTag = 10
    op.uniaxialMaterial('Elastic',flexMatTag,ut.Usmall)


    ParMatTag = 11
    op.uniaxialMaterial('Parallel',ParMatTag,matTag1,matTag2)

    #Define elements
    Transf1 = 1
    Transf2 = 2
    op.geomTransf('Linear', Transf1,1,0,0)  # Geometric transformation for rigid elastic elements 
    op.geomTransf('Linear', Transf2,0,0,1) 
    #Define rigid elastic elements
    Ae = 1000000000
    Ee =  210
    Ize = 8333333333300000
    Iye = 8333333333300000
    Ge = Ee/(2*(1+0.3))
    Je = 8333.333333300000

    eleTagiBeam = 1
    eleTagcBeam = 2

    op.element('elasticBeamColumn', eleTagiBeam, nodeTagCent, nodeTagi1, Ae, Ee, Ge, Je, Iye, Ize, Transf1)
    op.element('elasticBeamColumn', eleTagcBeam, nodeTagCent, nodeTagc1, Ae, Ee, Ge, Je, Iye, Ize, Transf1)

    #Doesn't work 
    #op.rigidLink('beam',nodeTagi1,nodeTagCent)
    #op.rigidLink('beam',nodeTagc1,nodeTagCent)

    eleTagis = 3
    eleTagcs = 4

    op.element('zeroLength', eleTagis, nodeTagi1,nodeTagi2,'-mat',ImpMatTag,ImpMatTag,rigidMatTag,rigidMatTag,rigidMatTag,rigidMatTag,'-dir',1,2,3,4,5,6)
    op.element('zeroLength', eleTagcs, nodeTagc1,nodeTagc2,'-mat',ConvMatTag,ConvMatTag,rigidMatTag,rigidMatTag,rigidMatTag,rigidMatTag,'-dir',1,2,3,4,5,6)


    eleTagSpoke = 100
    eleTagBP = 1000
    for i in range(0,nsp):
        op.element('elasticBeamColumn',eleTagSpoke+i,nodeTagCent,nodeTagCirc1+i, Ae, Ee, Ge, Je, Iye, Ize, Transf2)
        #op.rigidLink('beam',nodeTagCent,nodeTagCirc1+i)
        op.element('zeroLength',eleTagBP+i,nodeTagCirc2+i,nodeTagCirc1+i,'-mat',rigidMatTag,rigidMatTag,ParMatTag,flexMatTag,flexMatTag,flexMatTag,'-dir',1,2,3,4,5,6)


    #Run Gravity Analysis

    constTS = 1
    op.timeSeries('Constant', constTS)
    patternTagG = 1
    op.pattern('Plain',patternTagG,constTS)
    #op.load(nodeTagi2,0,0,0,0,0,0)
    #op.load(nodeTagc2,0,0,0,0,0,0)

    Tol = 1e-7 # convergence tolerance for test
    NstepGravity = 10
    DGravity = 1/NstepGravity
    op.integrator('LoadControl', DGravity) # determine the next time step for an analysis
    op.numberer('RCM') # renumber dof's to minimize band-width (optimization), if you want to
    op.system('BandGen') # how to store and solve the system of equations in the analysis
    op.constraints('Plain') # how it handles boundary conditions
    op.test('NormDispIncr', Tol, 100) # determine if convergence has been achieved at the end of an iteration step
    op.algorithm('Newton') # use Newton's solution algorithm: updates tangent stiffness at every iteration
    op.analysis('Static') # define type of analysis static or transient
    op.analyze(NstepGravity) # apply gravity
    op.loadConst('-time', 0.0) #maintain constant gravity loads and reset time to zero
    print('Static analysis completed')


    omega = []
    freq =  []
    T = []

    lamb = op.eigen('-fullGenLapack',4)

    for lam in lamb:
        omega.append((lam)**0.5)
        freq.append((lam)**0.5/(2*np.pi))
        T.append((2*np.pi)/(lam)**0.5)

    for j in  range(len(T)):
        print('T'+str(j+1)+' = '+str(T[j])+' s')


    #Define Recorder
    op.recorder('Node', '-file', direc+'/dispI_GM'+str(GM+1)+'.out','-time', '-node', nodeTagi2, '-dof', 1, 'disp')
    op.recorder('Node', '-file', direc+'/UpliftI_GM'+str(GM+1)+'.out','-time', '-node', nodeTagCirc1+4,'-dof',3,'disp')
    #op.recorder('Node', '-file', direc+'/reac.out','-time', '-node', nodeTagCirc2,'-dof','1','reaction')
    #op.recorder('Element', '-file', direc+'/force1.out','-time','-ele', eleTagBP, 'force')
    #op.recorder('Element', '-file', direc+'/force2.out','-time','-ele', eleTagBP+1, 'force')
    #op.recorder('Element', '-file', direc+'/force3.out','-time','-ele', eleTagBP+2, 'force')
    #op.recorder('Element', '-file', direc+'/force4.out','-time','-ele', eleTagBP+3, 'force')
    #op.recorder('Element', '-file', direc+'/force5.out','-time','-ele', eleTagBP+4, 'force')
    #op.recorder('Node', '-file', direc+'/disp5.out','-time','-node', nodeTagCirc1+4, '-dof', 3, 'disp')
    #op.recorder('Element', '-file', direc+'/force6.out','-time','-ele', eleTagBP+5, 'force')
    #op.recorder('Element', '-file', direc+'/force7.out','-time','-ele', eleTagBP+6, 'force')
    #op.recorder('Element', '-file', direc+'/force8.out','-time','-ele', eleTagBP+7, 'force')
    #op.recorder('Element', '-file', direc+'/forceI.out','-time','-ele', eleTagi, 'force')
    #op.recorder('Element', '-file', direc+'/forceC.out','-time','-ele', eleTagc, 'force')

    #vfo.createODB(model="Tank", loadcase="Pushover")

    #Run Dynamic Nonlinear time history analysis 
    print('Running NLTH Analysis, GM '+str(GM+1)+'...')
    GMfact = 1.0

    dta = dts/10
    TmaxAnalysis = dts*Npts
    GMfatt =  ut.g*GMfact

    IDLoadTagX = 100
    IDTimeSeriesX = 1000

    IDLoadTagY = 200
    IDTimeSeriesY = 2000

    #Define time series
    op.timeSeries('Path',IDTimeSeriesX,'-dt',dts,'-filePath',GMfile,'-factor',GMfatt)
    #op.timeSeries('Path',IDTimeSeriesY,'-dt',dts,'-filePath',GMfile,'-factor',GMfatt)

    #Define load pattern
    op.pattern('UniformExcitation',IDLoadTagX,1,'-accel',IDTimeSeriesX)
    #op.pattern('UniformExcitation',IDLoadTagX,2,'-accel',IDTimeSeriesY)

    Tol = 1e-8                          # convergence tolerance for test
    op.wipeAnalysis()
    op.integrator('Newmark', 0.5, 0.25) # determine the next time step for an analysis
    op.numberer('RCM')                  # renumber dof's to minimize band-width (optimization), if you want to
    op.system('BandGeneral')            # how to store and solve the system of equations in the analysis
    op.constraints('Plain')             # how it handles boundary conditions
    op.test('EnergyIncr', Tol, 50)      # determine if convergence has been achieved at the end of an iteration step
    op.algorithm('Newton')              # use Newton
    op.analysis('Transient') # define type of analysis static or transient

    ok =  op.analyze(10*Npts, dta)         # actually perform analysis; returns ok=0 if analysis was successful

    if(ok != 0):                    # analysis was not successful.
        # --------------------------------------------------------------------------------------------------
        # change some analysis parameters to achieve convergence
        # performance is slower inside this loop
        #    Time-controlled analysis
        ok = 0
        controlTime = op.getTime()
        while(controlTime < TmaxAnalysis and ok == 0):
            controlTime = op.getTime()
            ok = op.analyze(1, dta)
            if(ok != 0):
                print("Trying Newton with Initial Tangent ..")
                op.test('NormDispIncr', 1.0e-6, 100000, 0)
                op.algorithm('Newton', '-initial')
                ok = op.analyze(1, dta)
                op.test('EnergyIncr', Tol, 50)
                op.algorithm('Newton')
            if(ok != 0):
                print("Trying Broyden ..")
                op.test('NormDispIncr', 1.0e-6, 100000,  0)
                op.algorithm('Broyden', 8)
                ok = op.analyze(1, dta)
                op.algorithm('Newton')
            if(ok != 0):
                print("Trying NewtonWithLineSearch ..")
                op.algorithm('NewtonLineSearch', 0.8)
                ok = op.analyze(1, dta)
                op.algorithm('Newton')
            if(ok != 0):
                print("Trying KrylovNewton ..")
                op.algorithm(KrylovNewton) 
                ok = op.analyze(1, dta)
                op.algorithm('Newton')      

    print("Ground Motion Done. End Time:"+str(op.getTime()))
    op.wipe()


